home *** CD-ROM | disk | FTP | other *** search
/ Collection of Tools & Utilities / Collection of Tools and Utilities.iso / fortran / linpklib.zip / SPTSL.FOR < prev    next >
Text File  |  1984-01-06  |  3KB  |  92 lines
  1.       SUBROUTINE SPTSL(N,D,E,B)
  2.       INTEGER N
  3.       REAL D(1),E(1),B(1)
  4. C
  5. C     SPTSL GIVEN A POSITIVE DEFINITE TRIDIAGONAL MATRIX AND A RIGHT
  6. C     HAND SIDE WILL FIND THE SOLUTION.
  7. C
  8. C     ON ENTRY
  9. C
  10. C        N        INTEGER
  11. C                 IS THE ORDER OF THE TRIDIAGONAL MATRIX.
  12. C
  13. C        D        REAL(N)
  14. C                 IS THE DIAGONAL OF THE TRIDIAGONAL MATRIX.
  15. C                 ON OUTPUT D IS DESTROYED.
  16. C
  17. C        E        REAL(N)
  18. C                 IS THE OFFDIAGONAL OF THE TRIDIAGONAL MATRIX.
  19. C                 E(1) THROUGH E(N-1) SHOULD CONTAIN THE
  20. C                 OFFDIAGONAL.
  21. C
  22. C        B        REAL(N)
  23. C                 IS THE RIGHT HAND SIDE VECTOR.
  24. C
  25. C     ON RETURN
  26. C
  27. C        B        CONTAINS THE SOULTION.
  28. C
  29. C     LINPACK. THIS VERSION DATED 08/14/78 .
  30. C     JACK DONGARRA, ARGONNE NATIONAL LABORATORY.
  31. C
  32. C     NO EXTERNALS
  33. C     FORTRAN MOD
  34. C
  35. C     INTERNAL VARIABLES
  36. C
  37.       INTEGER K,KBM1,KE,KF,KP1,NM1,NM1D2
  38.       REAL T1,T2
  39. C
  40. C     CHECK FOR 1 X 1 CASE
  41. C
  42.       IF (N .NE. 1) GO TO 10
  43.          B(1) = B(1)/D(1)
  44.       GO TO 70
  45.    10 CONTINUE
  46.          NM1 = N - 1
  47.          NM1D2 = NM1/2
  48.          IF (N .EQ. 2) GO TO 30
  49.             KBM1 = N - 1
  50. C
  51. C           ZERO TOP HALF OF SUBDIAGONAL AND BOTTOM HALF OF
  52. C           SUPERDIAGONAL
  53. C
  54.             DO 20 K = 1, NM1D2
  55.                T1 = E(K)/D(K)
  56.                D(K+1) = D(K+1) - T1*E(K)
  57.                B(K+1) = B(K+1) - T1*B(K)
  58.                T2 = E(KBM1)/D(KBM1+1)
  59.                D(KBM1) = D(KBM1) - T2*E(KBM1)
  60.                B(KBM1) = B(KBM1) - T2*B(KBM1+1)
  61.                KBM1 = KBM1 - 1
  62.    20       CONTINUE
  63.    30    CONTINUE
  64.          KP1 = NM1D2 + 1
  65. C
  66. C        CLEAN UP FOR POSSIBLE 2 X 2 BLOCK AT CENTER
  67. C
  68.          IF (MOD(N,2) .NE. 0) GO TO 40
  69.             T1 = E(KP1)/D(KP1)
  70.             D(KP1+1) = D(KP1+1) - T1*E(KP1)
  71.             B(KP1+1) = B(KP1+1) - T1*B(KP1)
  72.             KP1 = KP1 + 1
  73.    40    CONTINUE
  74. C
  75. C        BACK SOLVE STARTING AT THE CENTER, GOING TOWARDS THE TOP
  76. C        AND BOTTOM
  77. C
  78.          B(KP1) = B(KP1)/D(KP1)
  79.          IF (N .EQ. 2) GO TO 60
  80.             K = KP1 - 1
  81.             KE = KP1 + NM1D2 - 1
  82.             DO 50 KF = KP1, KE
  83.                B(K) = (B(K) - E(K)*B(K+1))/D(K)
  84.                B(KF+1) = (B(KF+1) - E(KF)*B(KF))/D(KF+1)
  85.                K = K - 1
  86.    50       CONTINUE
  87.    60    CONTINUE
  88.          IF (MOD(N,2) .EQ. 0) B(1) = (B(1) - E(1)*B(2))/D(1)
  89.    70 CONTINUE
  90.       RETURN
  91.       END
  92.